eegmp_analysis/tools/[MEG]PLS/MEGPLS_PIPELINE_v2.02b/MEGPLS_COMPONENTS/MEGpipeline_NormMRI_Guts.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Normalises MRI and volumes via Fieldtrip and SPM. %
% Last modified: Jan. 15, 2014                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage:
%  [NormMRI, NormParams] = MEGpipeline_NormMRI_Guts...
%    (FTcfg, SPMcfg, MRIdata, Res, OutpathNormMRI)
%
% Inputs:
%  Note: FTcfg & SPMcfg can be generated & imported from a Builder .mat file.
%  FTcfg.WriteNormMRI  = Fieldtrip config for ft_volumewrite.
%  SPMcfg.Estimate     = SPM config estimate field from spm_get_defaults.
%  SPMcfg.Write        = SPM config write field from spm_get_defaults.
%
%  MRIdata        = Loaded FT structure containing MRI anatomy OR path to FT .mat MRIdata file.
%  Res            = Desired resolution (in mm) of output normalised MRI.
%  OutpathNormMRI = Output /path/filename of normalised MRI.


% Copyright (C) 2013-2014, Michael J. Cheung
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


function [NormMRI, NormParams] = MEGpipeline_NormMRI_Guts...
    (MainFTcfg, MainSPMcfg, MRIdata, Res, OutpathNormMRI)


NormMRI    = [];  % Initialize output variables
NormParams = [];


% Check number of input & output arguments:
if nargin ~= 5
    error('ERROR: Incorrect number of input arguments. See help for usage info.')
end
if nargout > 2
    error('ERROR: Incorrect number of output arguments. See help for usage info.')
end


% Load settings and check MRI:
SPMcfg.Estimate  = MainSPMcfg.Estimate;
SPMcfg.Write     = MainSPMcfg.Write;
SPMcfg.Write.vox = [Res, Res, Res];



%==========================%
% BEGIN MRI NORMALISATION: %
%==========================%

% If MRIdata is path to FT .mat, load it:
if ~isstruct(MRIdata)
    MRIdata = LoadFTmat(MRIdata, 'NormaliseSPM');
    if isempty(MRIdata)
        return;
    end
end


% Convert MRI into approximate RAS-orientation (SPM space):
% Keep initial transformation matrix (To emulate "ft_volumenormalise").
MRIdata = ft_checkdata(MRIdata, 'datatype', 'volume', 'feedback', 'yes');
MRIdata = ft_convert_units(MRIdata, 'mm');

orig	= MRIdata.transform;
MRIdata = ft_convert_coordsys(MRIdata, 'spm');
initial = MRIdata.transform / orig;


% Write volume fields to NIFTI:
[FolderNormMRI, NameNormMRI, ~] = fileparts(OutpathNormMRI);

disp('Writing MRI volume fields to NIFTI for SPM:')
VolumeFields = parameterselection('all', MRIdata);  % Get volume fields

for f = 1:length(VolumeFields)
    TempFilename{f} = [VolumeFields{f},'_',NameNormMRI];
    TempFilename{f}(strfind(TempFilename{f}, '.')) = [];  % ft_volumewrite cannot read dots.
    
    if strcmp(VolumeFields{f}, 'anatomy')
        TempAnatFile = [TempFilename{f},'.nii'];
    end
    
    FTcfg.WriteNormMRI			 = [];
    FTcfg.WriteNormMRI			 = MainFTcfg.WriteNormMRI;
    FTcfg.WriteNormMRI.parameter = VolumeFields{f};
    FTcfg.WriteNormMRI.filename  = TempFilename{f};
    FTcfg.WriteNormMRI.filetype  = 'nifti';
    
    ft_volumewrite(FTcfg.WriteNormMRI, MRIdata)
end

MRIdata = [];  % Free memory


% Normalise anatomy volume and get parameters:
Template = [spm('dir'),'/templates/T1.nii'];

disp('Normalising anatomy:')
NormParams = spm_normalise(Template, TempAnatFile, [], [], [], SPMcfg.Estimate);
spm_write_sn(TempAnatFile, NormParams, SPMcfg.Write);


% To emulate "ft_volumenormalise:" 
% Determine the affine source -> template coordinate transformation:
TemplateHdr = spm_vol(Template);
AnatHdr     = spm_vol(TempAnatFile);
Final = TemplateHdr.mat * inv(NormParams.Affine) * inv(AnatHdr.mat) * initial;

TemplateHdr = [];  % Free memory
AnatHdr     = [];


% Read normalised anatomy back into FT:
NormMRI               = ft_read_mri(['w',TempAnatFile]);
NormMRI.params        = NormParams;
NormMRI.initial       = initial;
NormMRI.coordsys      = 'spm';
NormMRI.unit          = 'mm';
NormMRI.cfg.final     = Final;
NormMRI.cfg.spmparams = NormParams;

system(['rm ',TempAnatFile]);  % Remove remnant volume field files.
system(['rm w',TempAnatFile]);
Final = [];  % Free memory


% Normalise other volume fields and read back into FT:
disp('Apply normalisation to other fields:')
for f = 1:length(VolumeFields)
    if strcmp(VolumeFields{f}, 'anatomy')
        continue;  % anatomy already done
    end
    
    TempFile = [TempFilename{f},'.nii'];
    spm_write_sn(TempFile, NormParams, SPMcfg.Write);
    
    hdr     = spm_vol_nifti(['w',TempFile]);
    img     = spm_read_vols(hdr);
    NormMRI = setsubfield(NormMRI, VolumeFields{f}, img);
    
    system(['rm ',TempFile]);  % Remove remnant volume field files.
    system(['rm w',TempFile]);
    hdr = [];  % Free memory
    img = [];
end

if isfield(NormMRI, 'inside')  % Convert inside back to logical
    NormMRI.inside = abs(NormMRI.inside-1)<=10*eps;
end


% Align volume (voxel-space) with headcoordinate axes (SPM world-space):
% Note: Yields identical results as loading & saving with NIFTImatlab.
NormMRI = align_ijk2xyz(NormMRI);


% Save NormMRI:
CheckSavePath(OutpathNormMRI, 'NormaliseSPM');
save(OutpathNormMRI, 'NormMRI');